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Abstract. We apply ftp-cloud method to the radial Dirac eigenvalue problem. The difficulty of 
occurrence of spurious eigenvalues among the genuine ones in the computation is resolved. The 
method of treatment is based on assuming ftp-cloud Petrov-Galerkin scheme to construct the 
weak formulation of the problem which adds a consistent diffusivity to the variational formula- 
tion. The size of the artificially added diffusion term is controlled by a stability parameter (r). 
The derivation of r assumes the limit behavior of the eigenvalues at infinity. The parameter r is 
applicable for generic basis functions. This is combined with the choice of appropriate intrinsic 
enrichments in the construction of the cloud shape functions. 



1. Introduction. 

In the last decades, several numerical methods have been derived to compute the eigenvalues 
of operators. The need of accurate computation of eigenvalues is considered due to their signifi- 
cant applications in many disciplines of science. Mathematically, if a matrix or a linear operator 
is diagonalized, then by the spectral theorem, it can be analyzed by studying its corresponding 
eigenvalues, i.e., transforming the matrix or operator to a set of eigenfunctions which can be eas- 
ily studied. From the physical point of view, the eigenvalues possess a wide range of information 
about the behavior of the system governed by an operator. This information might be all what 
is needed to answer many questions regarding the system properties such as stability, positivity, 
boundedness, asymptotic behavior, etc. In mechanics, eigenvalues play a central role in several 
aspects such as determining whether the automobile is noisy, whether a bridge will collapse by 
the water waves, etc. Also, the eigenvalues describe how the quantum state of a physical system 
changes in time (Schrodinger equation). They also represent the electrons relativistic energies 
and describe their motion in the atomic levels, this is the well-known Dirac equation, which is 
the core of the present work. 

The calculation of energy levels in Helium-like ions, where the interaction between two elec- 
trons takes place, can be determined by studying the electrons correlation which is part of quan- 
tum electrodynamic effects (QED-effects). A scheme for calculating QED-effects p^EHlGElllQ] is 
based on a basis set of relativistic Hydrogen-like ion wave eigenfunctions (of the Dirac operator) . 
Meanwhile, the numerical computation of the Dirac operator eigenvalues encounters unphysical 
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values (do not match the physical observations) called spurious eigenvalues or spectrum pollu- 
tion. The spurious eigenvalues result in rapid oscillations in the wave functions, hence, in many 
cases, affecting the computation reliability of the basis set in the practical atomic calculations. 

The spurious eigenvalues are an effect of the numerical methods and are found in the com- 
putation of many problems other than the Dirac eigenvalue problem [TJ [2l [37] SI] • For general 
eigenvalue problems, spurious eigenvalues are reported in [36]. The occurrence of the spuriosity 
is related to mismatching of desired properties of the original solution in the numerical formu- 
lation. Also in the computation of electromagnetic problems the spuriosity is perceived [344 139], 
Two leading approaches are derived to solve this difficulty; Shabaev et al. [UJ have related 
the spuriosity to the symmetric treatment of the large and small components of the Dirac wave 
function. Their approach, for removing the spurious eigenvalues, is based on deriving dual 
kinetic-balance (DKB) basis functions for the large and small components. Almanasreh et al. 
[2] have allied the occurrence of spurious eigenvalues to the incorrect treatment of the trial and 
test functions in the finite element method (FEM). They proposed a stability scheme based on 
creating diffusivity by modifying the test function so that it includes a gradient-based correction 
term. 

In this work, we apply /ip-cloud method [151 I37j to the radial Dirac eigenvalue problem. The 
technique is known as one of the meshfree methods (MMs) [U [181 130] [31] [35] that allows two 
different enrichments, intrinsic and extrinsic, to be built in the construction of the basis functions. 
The method was previously applied for several problems, e.g., the Schrodinger equation |10j . 
Mindlin's thick plate model [19], and Timoshenko beam problems [32], etc. Here, we apply 
/ip-cloud method based on the Galerkin formulation. This means that it is required to evaluate 
the integrals in the weak formulation of the particular equation, thus a background mesh must 
be employed in the integration techniques. Therefore, the /ip-cloud method used here is not 
really a truly MM. However, all other features of MMs are maintained in our approximation. 

In order to treat the spuriosity problem, we stabilize the computation by considering instead 
an /ip-cloud Petrov-Galerkin (/ip-CPG) method which may be considered as a new technique of 
the general meshfree local Petrov-Galerkin (MLPG) [3[ 1171 [28] . The stability scheme is based 
on adding consistent diffusion terms without changing the structure of the equation. The size 
of the additional diffusivity is controlled by a stability parameter. 

Consider the radial Dirac eigenvalue problem H K §{x) = A<£(x), where <5(x) = (F(x), G(x)) t 
is the radial wave function, A is the electron relativistic energies (eigenvalues), and H K is the 
radial Dirac operator given by 

/ mc 2 + V(x) c(-D x + %) \ 
Hk \ c(D x + 1) -mc 2 + V(x) J ■ 

The constant c is the speed of light, m is the electron mass, V is the Coulomb potential, 
D x = d/dx, and k is the spin-orbit coupling parameter defined as k = (— l) 3+e - + 2 (j + i), where j 
and I are the total and the orbital angular momentum quantum numbers respectively. The weak 
formulation of the problem is to find A G M and $ in a specific function space such that for every 
test function in some suitable space we have J Q ^ l H K ^dx = A j n ^Qdx. The usual /ip-cloud 
Galerkin approximation is to let to be {ip, 0)* and (0, ip) 1 , where tp is in the same space as of the 
two components of $. To discretize the weak form, the components of the trial function $ and 
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the test function tp are chosen from a finite set of /ip-cloud basis functions which are constructed 
using moving least-squares method. Since the radial Dirac operator is dominated by advection 
(convection) terms, the /ip-cloud approximation will be upset by spurious eigenvalues. 

To stabilize the /ip-cloud approximation, the hp-CP G method is used. In this formulation, the 
test function \& is assumed to belong to a function space different from that of the trial function 

in the sense that is chosen in the form (ip, rip') and (tiJj' , ip) 1 where ip belongs to the same 
space as the two components of <J>. The correction term rip' is used to add artificial viscosity, 
controlled by r, to stabilize the convection terms. The derivation of the stability parameter r 
follows the principle used in [2] for the FEM. Two assumptions are considered in deriving r; (i) 
the operator limit as the radial variable x tends to infinity, thus obtaining an approximation of 
the limit point of the eigenvalues (depending on r) which can be compared to the theoretical 
limit point eigenvalue [20] . (u) considering the dominant terms with respect to the speed of light 
(c). 

The paper is organized as follows; in Section 2, some preliminaries about the Dirac equation 
are presented, also we shed some light over the occurrence of the spuriosity. In Section 3, 
the construction of the /ip-cloud functions is provided, also coupling with the FEM to impose 
essential boundary conditions (EBCs) is explained. The hp-CPG method and the derivation of 
the stability parameter are treated in Section 4. In the last section, Section 5, we present some 
numerical results and provide necessary discussion about the stability scheme. 

2. The radial Dirac eigenvalue problem and the spuriosity 
The free Dirac space-time equation is given by 

d 

(1) ih— u(x,t) =H u(x,t), u(x,0) = u (x), 

where K is the Planck constant divided by 2vr, and H : if 1 (R 3 ; C 4 ) — > L 2 (M 3 ;C 4 ) is the free 
Dirac operator acting on the four-component vector u, given by 

(2) H = -i hca ■ V + mc 2 f3 . 
The 4x4 Dirac matrices a = (ai,a2) a 3) an d /3 are given by 

Here / and are the 2x2 unity and zeros matrices respectively, and <7j's are the 2x2 Pauli 
matrices 

ai= {°i J)' a2= o 1 )' anda3= (i -i)- 

Note that separation of variable yields the Dirac eigenvalue problem, i.e., by assuming u(x, t) = 
u(x)6(t) in (P) one gets 



(3) 



Hqu(x) = Au(x). 
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The operator Hg is self adjoint on i7 1 (]R 3 ; C 4 ), it describes the motion of the electron that moves 
freely without external forces. The free Dirac operator with Coulomb potential is given as 

(4) H = H + y(x)7, 

where V(x) = t^-, and I is the 4x4 unity matrix. The spectrum, denoted by a, of the 
Coulomb-Dirac operator is <r(H) = (— oo,— mc 2 ] U {Afc}fe £ N U [mc 2 , +oo), where {Afc}fc £ N is a 
discrete sequence of eigenvalues in the gap (— mc 2 , mc 2 ) of the continuous spectrum. 

The radial Coulomb-Dirac operator (radial Dirac operator) can be obtained by separation of 

1 / F(x)x K ,m(TZ>, 6) 



variables of the radial and angular parts, i.e., by assuming urx) - 

where x represents the radial variable. The spherical symmetry property of the angular function 
X is the key point in the derivation of the radial part. Let 3>(cc) = (F(x), G{x)) t , the radial 
Dirac eigenvalue problem is given as 

(5) H K $(x) = \$(x), where 

[ j V c(D x + f) -mc 2 + V{x) 

The radial functions F{x) and G{x) are called respectively the large and small components of 
the wave function <&{x). 

The well-known difficulty of solving the radial Dirac eigenvalue problem numerically is the 
presence of spurious eigenvalues among the genuine ones that disturb the solution and conse- 
quently affect the reliability of the approximated eigenstates. Here we follow [2] for the clas- 
sification of the spurious eigenvalues; the first category is the so-called instilled spuriosity, and 
the second category is the unphysical coincidence phenomenon. The first type is those spurious 
eigenvalues that may occur within the true eigenvalues (i.e., they occur between the genuine en- 
ergy levels), this type of spuriosity appears for all values of k. The second type is the unphysical 
assigning of almost the same first eigenvalue for 2s 1 / 2 (k = —1) and 2pi/ 2 {n = 1), 3p 3 / 2 (fc = —2) 
and 3g?3/2(k = 2), 4c£5/ 2 (k = —3) and 4/5/2 (k = 3), and so on. This phenomenon is rigorously 
studied in [32] from theoretical aspect and in [33] from numerical viewpoints. 

Apparently, most authors [TJ [2j [371 ETJ agree that the incorrect balancing and the symmetric 
treatment of the large and small components of the wave function or of the test and trial 
functions in the numerical methods are the core of the problem. In the present work, we relate 
the occurrence of spuriosity, of both categories, to the unsuitable function spaces and to the 
symmetric treatment of the trial and test functions in the weak formulation of the equation. To 
clarify, we rewrite ([5]) to obtain an explicit formula for each of the two radial functions, so by 
defining w ± {x) = ±rac 2 + V(x) we have, see [2], 

(7) F"(x) + ll (x,\)F'(x)+ l2 {x,\)F{x) = 0, 

(8) G"{x) + 0i (x, X)G'(x) + 9 2 (x, X)G(x) = , 
where 

7i(x,A) = —— -, 6>i(x,A)- 



w (x) — A ' w + (x) — A ' 
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72 (as, A) 

and 



(w + (x) - A) (w~ (x) - A) k 2 + k kV'{x) 



x 2 x(w~(x) — A) 



, XN (w + (x) - A) (w~(x) - A) k 2 -k 
jfo A) = ^ ^ q ^ >- =- + 



c 2 x 2 x(w + (x) — A) 

According to (J7|) and ([8]) , the components F and G should be continuous and have continuous 
first derivatives. Thus, the suitable choices of function spaces for the computation of the radial 
Dirac eigenvalue problem are those that possessing these properties. Assuming appropriate 
spaces helps in partial elimination of spurious eigenvalues, and does not help overcoming the 
unphysical coincidence phenomenon. In [2J, the same argument is accounted, where the FEM is 
applied to approximate the eigenvalues using linear basis functions. 

Table 1. This table, taken from [2], shows the first computed eigenvalues of the 
electron in the Hydrogen atom. 



Level 


K = 1 


K = -1 


Rel. Form, k = — 1 


1 


-0.50000665661 


-0.50000665659 


-0.50000665659 


2 


-0.12500208841 


-0.12500208839 


-0.12500208018 


3 


-0.05555631532 


-0.05555631532 


-0.05555629517 


^ 


-0.03141172061 


-0.03141172060 


Spurious Eigenvalue 


4 


-0.03118772526 


-0.03118772524 


-0.03125033803 


5 


-0.01974434510 


-0.01974434508 


-0.02000018105 



In Table 1, 400 nodal points are used to discretize the domain, and the computation is 
performed for point nucleus. The shaded value in the first level is what meant by the unphysical 
coincidence phenomenon, and the two shaded values after the third level are the so-called instilled 
spuriosity. If the basis functions are chosen to be C 1 -functions, then some instilled spurious 
eigenvalues are avoided as indicated in [2]. Therefore, after applying the boundary conditions, 
homogeneous Dirichlet boundary condition is assumed for both radial functions, the proposed 
space is JC(f2) := C 1 (0) n Hq(Q). Also, the radial functions are mostly like to vanish at 
the boundaries in a damping way (except some states), consequently homogeneous Neumann 
boundary condition should be taken into account. The exceptional states are ls]/2 and 2pi/2, 
in this case at the upper boundary the same treatment is considered as of the other states, but 
the first derivative of these two states at the lower boundary is not zero. Here we will assume 
general boundary condition, that is, homogeneous Dirichlet boundary condition. Thus, from 
now on, the space !K(f2) is considered. 

What deserves to dwell upon is that in numerical methods applied to convection dominated 
problem, the solution is disturbed by spurious oscillations. This instability gets worse as the 
convection size increases. The following two numbers are considered as tools to measure the 
dominance of the convection term 



(9) 
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where Pej and Daj are known as the grid Peclet and Damkohler numbers respectively, hj is 
the size of the element interval Ij , Uj and sj are respectively the coefficients of the convection 
and the reaction terms corresponding to and K is the diffusivity size. The difficulty arises 
when either the convection coefficient or the source term is larger than the diffusion coefficient, 



convection dominated one. 

For the simplified equations ([7]) and (JSj) -> it is easy to check that 2PeDa admits very large 
values if small number of nodal points is considered regardless the sizes of |A|, Z, and k. Even 
with mesh refinement, 2PeDa still admits very large values at some positions (2PejDaj >> f 
for some j). The Peclet number, Pe, for the equation that involves the function F, is always 
less than 1. But for the equation that corresponds to G, Pe admits a value greater than one, 
exactly at the nodal point where Uj changes its sign, here refining the mesh does not bring Pe 
below one for all choices of A, Z, and k. Hence, ([7]) and (|8j) are dominated by convection terms. 
Thus the approximated solutions F and G, will be upset by unphysical oscillations. The draw 
back is that these oscillations in the eigenfunctions result in some unphysical eigenvalues, the 
spurious eigenvalues. For detailed study on convection dominated problems we refer to [71 [8]. 



To build the /ip-cloud functions, MLS method is applied which allows easily p-enrichment to 
be implemented and to desired fundamental characters to be enriched in the approximation. 
MLS is a well-known approximation technique for constructing meshfree shape functions in 
general. It applies certain least square approach to get the best local approximation, then the 
local variable is let to 'move' to cover the whole domain. 

Consider an open bounded domain Q C M with boundary dQ, assume X = {x\,X2, . . . ,x n } is 
a set of n arbitrary points in 0. Let W = {wi}^ =1 be a set of open covering of defined by X 
such that Wi is centered at Xi and C U™ =1 Wj. 

Definition 1. A set of functions {V , i}£ = i is called a partition of unity (PU) subordinated to the 
cover W if 

(1) Yli=i iM^O = 1 f or a M x ^ ^- 

(2) tp E Cq^Wi) for i = 1, 2, . . . , n, where s > 0. 

Let P = {pi(x),p2(x), . . . ,Pm(x)} be a set of basis of polynomials (or any basis of suitable 
intrinsic enrichments). Suppose that is a continuous function defined on Q and that its 

values, ^i, at the points Xi G O, i = 1, 2, . . . , n, are given. To approximate ty(x) globally by 
^h(x), firstly fy(x) is approximated locally at x € by defined in terms of the given basis 
set P as 



i.e., when Pej > 1 or when 2PejDaj 




> 1 respectively, then the associated equation is a 



3. Moving least-squares (MLS) approximation 



(10) 



Js&ix) = P\x)a{x) 



where t denotes the usual transpose. The unknown coefficients a(x) are chosen so that Jx*& is 
the best approximation of \t in a certain least squares sense, this is achieved if a is selected to 
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minimize the following weighted least squares L 2 -error 



[11) 



4(a) = V^(^— ^i)(P*(xi)a(x) - * 4 ) 2 



t=l 



ft 



where is a weight function introduced to insure the locality of the approximation, and pi 
is the dilation parameter which controls the support radius of ifi at x%. To minimize I x with 
respect to the vector a, the first derivative test is applied, i.e., we set Q^. = q which gives 



dl x . 
7^7 = 2^ 



da; 

J i=l 



Pi 



)2p j (x i ){P t {x l )a(x)-^i) = 0, j = 1,2,..., m. 



The above system can be written as 

(12) M{x)a{x) - B(x)V = 0, 

where M{x) = 3>'T(x):P, B{x) = !P*T(x), = [*i,¥ 2 , ■ ■ ■ , and a*(x) = [a^x) , a 2 (x) , . . . 
a m (x)], with J 5 and T(x) defined as 



/ Pi(xi) p 2 (xi) ... p TO (xi) \ 
Pi(x 2 ) p 2 (x 2 ) ■■■ Pm{x 2 ) 

\ Pl(x n ) p 2 (x n ) ... p m {Xn) ) 



and 



T(x) 



V 



o 












We proceed from equation (|T2|) to get 

(13) a(x) = M _1 (rc)i?(x)^. 

Thus 

J £ *(x) = P\x)a(x) = P*(x)M~ 1 (x)S(x)^. 

The global approximations is then obtained by assuming x arbitrary, i.e., by letting x move over 
the domain, viz, the solution is globalized by considering ^(x) ~ lim J x ^(x) =: ^h(x), thus 



(14) 



i=l 



with ipi(x) = P t (x)M 1 (x)Bi(x), and Bi(x) = (fi(^-^ :L )P(xi). To sum up, ^ can be written as 



(15) 



n 1 n 

= p\x) ( V Vl (^i)p( Xt )p\ Xl )) ~ V ^^—^)P{ Xl )^ t . 



i=i 



pi 



The first derivative of fa is given by fa jX = = P t x M- x B i -P t M- x M x M- x B i ^P t M- 1 B i ^. 

Below we shall need the consistency concept. 
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Definition 2. A set of functions {iii(x)} is consistent of order m, if Ui(x)ty(xi) = *$(x) for 

i 

all where *$(x) = \$\ < m}. 

To increase the order of consistency of the approximation, the complete representation of the 
/ip-cloud functions consists of the set of PU functions tpi (x) and monomial extrinsic enrichment 
basis functions P as 

n no n no 

i=l j=i t=i j=i 

Note that P can be any type of basis functions, but the most used is monomials since they 
provide good approximation for smooth functions. The monomials Pj(x), according to [47] . 
should be normalized by the measure of the grid size at Xj to prevent numerical instability. 
Nevertheless, in applying the /ip-cloud approximation for the radial Dirac eigenvalue problem, 
we will use a stability scheme based on the MLPG method, for that we will not be interested 
in concerning extrinsic enrichments in the computation (P = {1}, a monomial of degree zero). 
The point of this setting follows [I], where six different realizations of MLPG restricted only to 
intrinsic enrichment basis are considered. It is found that extrinsic enrichments in the MLPG 
method cause numerical stability problems, because the behavior of their derivatives has large 
oscillations, which is not the case in the usual MMs. Hence, in the present work, only intrinsic en- 
richments, P{x), are considered, and thus the approximation with the /ip-clouds is given by (p2 



The weight function <p>i plays the most important role in the definition of the /ip-cloud shape 
function, it is defined locally on the cover wi around X\. The function ipi can also be chosen the 
same for all nodes, in this case we write cpi = <p, which is the case we consider in this work. The 
hp-cloud, tpi, inherits the properties of the weight function ip such as continuity and smoothness. 
In other words, if ip is continuous with continuous first derivative, then so is tpi, provided that 
the continuity of the enrichment basis P(x) and its first derivatives is ensured. As for the Dirac 
large and small components, F and G, the proposed space is "K, therefore, the weight function 
ip should be at least C -function. For this purpose, we will consider quartic spline (which is a 
C 2 -function, sufficiently enough) as a weight function defined by 

1 -6r 2 + 8r 3 - 3r 4 , r < 1, 
0, r>l, 



(16) <p(r) 



where r = ^ x *' . 

Pi 

The set functions {ipi}f =1 builds a PU, also the set of their first derivatives {ipi, x }i=i builds 
a partition of nullity (PN) (^2i = i' l Pi,x( x ) = for all x £ Q), see Figure 1. The computational 
effort of evaluating the integrals in the weak form of the /ip-cloud approximation is more time 
consuming compared to mesh-based methods (the shape functions are of the form <pi only), this 
is due to the fact that the derivative of the shape function tpi tends to have non-polynomial 
characters, also due to the time needed for matrix inversion in evaluating the shape functions. 

Since the Kronecker delta property being not a character of tpi (tpi(xj) ^ 8ij), then at each 
node there are at least two nonzero shape functions. Thus, to have the value of the approximated 
function at a node, all nonzero shape functions values at that node should be added. The missing 
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Figure 1. PU %>-clouds (to the left) and their PN first derivatives (to the right). 
Quartic spline is used as a weight function. 



of the Kronecker delta property causes a problem in imposing EBCs, and thus other techniques 
are used to solve this difficulty, see below. 

The intrinsic enrichment P(x) has an important effect in the definition of the /ip-cloud func- 
tions. All known fundamental characters, such as discontinuities and singularities, about the 
sought solution can be loaded on the intrinsic functions. Consequently, more time is saved; it 
is not needed, in general, to assume very large number of nodal points to capture a desired 
behavior of the approximated function while most of the solution features are inserted in the 
approximation itself. On the other hand, stability is enhanced particularly when there are some 
crucial characters that can not be captured by usual numerical methods, for example solving 
equations with rough coefficients that appear, e.g., in composites and materials with micro- 
structure, problems with high oscillatory solutions, or eigenvalue problems that admit spurious 
solutions in the computation of the discrete spectrum. 



Imposition of essential boundary conditions (EBCs) 



The radial Dirac eigenvalue problem assumes homogeneous Dirichlet boundary condition, 
while it is known that the /ip-cloud approximation (MMs in general) can not treat this condition 
naturally, this is because the lack of the Kronecker delta property of the shape function. This 
is in contrast with most mesh-based methods, where the basis functions admit this property, 
and thus applying EBCs is straightforward (as in FEM) by omitting the first and the last basis 
functions. 

In MMs in general, the widely applied techniques for imposing EBCs are Lagrangian multipli- 
ers, penalty condition, and coupling with finite element shape functions. Lagrangian multiplier 
is a very common and accurate approach for the imposition of EBCs. The disadvantage of this 
technique, see e.g. [Ml [45], is that the resulted discrete equations for a self-adjoint operator are 
not positive definite (contains zero at the main diagonal) nor banded. Also the structure of the 
system becomes awkward, i.e., instead of having M as a resulting matrix of discretization of the 

/ M Im \ 

Galerkin formulation, the system [ is obtained, where Im is the EBC-enf or cement 

I Im I 
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vector. EBCs can also be imposed by penalty condition |18^ I36j . the problem of applying this 
technique is the negative effect on the condition number of the resulting discrete equations. 

The most powerful and safe method to enforce EBCs is coupling MMs with the FEM, applied 
for the first time in [26]. With this approach, the meshfree shape functions of the nodes along 
boundaries are replaced by finite element basis functions. In one dimensional case, the fop-cloud 
functions at the first two and the last two nodes are replaced by finite element functions, and 
thus EBCs are, as in the FEM, simply imposed through eliminating the first and last added 
finite element functions, see Figure 2. 



hp-cloud and finite element coupled functions for exponentially distributed nodal points hp-cloud and finite element coupled functions for exponentially distributed nodal points 




domain O domain £1 



Figure 2. Coupled fop-cloud and finite element functions: general coupling (to 
the left), and coupling for the purpose of imposing EBCs (to the right) (two 
finite element functions are sufficient). Linear (hat) functions are used as finite 
element functions, and quartic spline as a weight function in the fop-clouds. 

Two efficient approaches of coupling MMs with the FEM are coupling with Ramp functions 
[5], and coupling with reproducing conditions [23]. With the former approach, the derivative of 
the resulting coupled approximation function at the boundary of the interface region, 17 tsn in 
Figure 2, is discontinuous and the consistency is of first order. To ensure the continuity of the 
first derivative of the coupled function and to obtain consistency of any order, we consider the 
latter approach. Using MLS approximation method as before, the approximation resulting from 
coupling fop-cloud and finite element functions with the reproducing conditions is given as (see 

e.g. m) 

where Si are the finite element shape functions, and M is as defined before. From Figure 2, 
it can be seen that finite element functions are only complete in fi FEM , and that in il MM only 
fop-clouds are present. In the transition interface region, Q tsn , the existence of incomplete finite 
element functions modifies the existed fop-clouds there, and thus coupled fop-cloud and finite 
element functions are obtained. 
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4. The scheme and the stability parameter 



Since the radial Dirac eigenvalue problem is convection dominated, the /ip-cloud approxima- 
tion of it will be unstable. As most of applications of numerical methods, certain modifications 
are used to stabilize solutions [21 [31 [71 [H [121 [25] • Therefore, instead of considering the fop-cloud 
approximation for the radial Dirac eigenvalue problem, we will apply the hp-CPG method to cre- 
ate diffusion terms to stabilize the approximation. The hp-CPG method is a consistent method 
in the sense that the solution of the original problem is also a solution to the weak form. The 
size of the added diffusivity is controlled by a stability parameter. To set the scheme, consider 
the radial Dirac eigenvalue problem H K Q = A<£, the usual fop-cloud method is formulated by 
multiplying the equation by a test function and integrating over the domain f2 



[17) 



WH^dx = A / W^dx. 
n Jn 



To discretize (fTTj) let ^ be defined as (tpi, 0)* and (0, tpiY, i = 1, 2, . . . , n, where tpi is the fop-cloud 
basis function and 



F(x) 
G(x) 



(18) $(x) = 

The elements fj and gj are the nodal values of F and G respectively. This yields 



Y,]=i9j^j( x ) 



(19) ^(w + (x)V>j(x) , ^iix^fj+^i-cip'jix) + —tpj{x) , ij)i{x))gj = A^(V>j(z) , ipi(x))fj 

3=1 j=l J'=l 

and 



j=i j=i j=i 



the bracket (• , •) is the usual L 2 (Q) scalar product. After simplifying, equations ([T9|) and ([20 
lead to the symmetric generalized eigenvalue problem 

(21) AX = XBX . 

The block matrices A and B are defined by 



mc 2 M 00 o + 


-cMoio + ckM 001 


cMqio + creMooi 


-mc 2 M 000 + M ^ 



(22) A 



where M^ st are n x n matrices given as 



and B 














(23) 



(m; 



rstJi] 



^ S) 4 r) x- t q(x)dx , (>)(x) = ^{x)) • 



The vector X is the unknowns defined as (/ , g) , where / = (/i, /2, • • • , fn) and g = (gi, g2, ■ ■ ■ , 9n)- 

To formulate the fop-CPG method, the test function \& is modified to include the first derivative 
of the basis function in order to introduce the required diffusivity. This leads to assume \E r as 
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(tpjTtp'Y and (rip' , ip) 1 in ()17|) . where r is the stability parameter that controls the size of the 
diffusion terms, ip = ipi, and the functions F and G are given by (j!8p . thus we get 

(24) (w+F , V) + i-cG' + ^G,iP) + (Re 2 (F, G) , t^j 1 ) = X(F , ^) 
and 

(25) (cF' + —F, tp) + {w~G , tp) + (iie 1 (F, G) , r^') = A(G , . 
The residuals Fe 1 (F, G) (x) and i?e 2 (F, G) (x) are defined as 

(26) Re 1 (F, G) (x) = {W + F - cG' + —G) (x) , 



(27) Re 2 (F, G) (x) = (?"G + cF' + — F) (x) , 

where W ± (x) = w^(x) — A. This results in the usual /ip-cloud approximation with addition to 
perturbations sized by r as follows 

(28) AX = XBX . 

The perturbed block matrices, A and B, are respectively in the forms A = A + tA and B = 
B + t23, where A and B are given by 



(29) A 



cMno + ckMioi 


-mc 2 M 100 + M^ ' 


mc 2 M 100 + M^ 


-cMno + ckMiqi , 



and B 






Mioo 


Mioo 






The system f)28f) is not symmetric, thus complex eigenvalues may appear if the size of r is 
relatively large. In the FEM, an explicit representation for r is obtained [2], where the basis 
functions have the Kronecker delta property, hence the basis functions have regular distribution 
along the domain and only the adjacent basis functions intersect in one and only one subinterval. 
Thus the resulted system consists of tridiagonal matrices, this makes the derivation of r easier 
and an explicit representation is feasible. In MMs in general, a basis function is represented 
by cloud over a nodal point, with domain of influence, p, that may cover many other nodal 
points. So the resulting matrices can be filled with many nonzero elements, hence the number 
of non-vanishing diagonals in these matrices is arbitrary (greater than 3) and depending on the 
size of p. Therefore, we can not write an explicit representation for r that depends only on a 
given mesh. Instead, r will be mainly represented by some of the computed matrices obtained 
from the usual /ip-cloud method. 

The derivation of r assumes the limit Dirac operator in the vicinity of x at infinity. This 
presumable simplification is inevitable and justifiable; the derivation leads to an approximation 
of the limit point eigenvalue which depends on r, where, in [20], the theoretical limit is proved 
to be mc 2 , hence we can minimize the error between the theoretical and the approximated limits 
to get r. By considering the limit operator at infinity, we consider the troublesome part (that 
includes the convection terms) of the operator which is mostly needed to be stabilized. Besides 
that, one is obliged to assume that the stability parameter should be applicable at all radial 
positions x € f2, particularly the large values of x. 
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Theorem 1. Let Mqoo & n d -^100 be the n x n computed matrices given by \23\i. and let o~ji and 
r\ji be the corresponding entries respectively. Define $ as 



(30) 



-J2 hk > 


i < j 


k=i+l 




o, 


i = j 


i 


i > j 



where hk is the displacement between the nodes x k and x k _\. Then the stability parameter, Tj, 
for an arbitrary j th row of the matrices in A and 23 is given by 



(31) 



Tj — Yj a ji$ji / Y2 ty^ji 



i=l i=l 

Proof . Consider the limit operator of the radial Dirac eigenvalue problem in the vicinity of x 



at infinity 
(32) 



mc 2 —cD x 
cD x —mc 2 



F(x) 
G(x) 



X 



F(x) 
G(x) 



The hp-CPG variational formulation of (|32j) (which is equivalent to assume a limit passage as 
x — > oo of the equations (f24"|) and ([25]) ) provides 



(33) 
and 
(34) 



(mc 2 - A)M 000 / + TcM 110 f - (rmc 2 - c + t\)M 100 9 = 
(rmc 2 - c - rA)Mioo/ - TcM uo g - (mc 2 + X)M 000 g = , 



where, as defined before, / = (/i,/ 2 , ■--,/«) and g = (gi,g 2 , ■ ■ -,g n )- Let a k , r) k , and g k , for 
k = 1, 2, . . . , n, be the corresponding j th row entries of -Mootb -^100) and Muq respectively. To 
obtain Tj, we consider the j th rows in (|33p and (|34p . this together with the Lemma 1 below gives 



(35) 



it it, n it 

(mc 2 -X^j ( Y Vkfj+Y, ak ( ma&k + ^ k / c ^) 9 i) +TC ( Y Skfj + Y Sk(mc& k + 

k=l k=l k=l k=l 

n n 

+(fi k /c)\)gj) - (rmc 2 -c + tXJ ( + ^Vk (mcd k ~ ('&k/c)X)f j \ = 



k=l 



k=l 



and 
(36) 



rmc — c 



rXj( Y^kfj + Y2 llk ( mc ® k + (^ k / c ) X )9j) ~ T c (Y 6k9 ^ + 



+ Yj P k ( mC ^ k - ($k/c)X)fj 



k=l 



(mc 



k=l 

n 



k=l 



+ ^){Y ak9 j +X) crfe ( mc1?fe ~ ^k/c)X)fj \ = 



k=l 



k=l 



k=l 
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Lemma 1. Let f i and g; t be respectively the i th nodal values of F and G of the limit equation 
(CHJj. Freeze j, and let $j be given by \30\) for the given j. Then for i = 1, 2, . . . , n 



finf j +(mc& i + (0 i /c)\)g j . 
gj + (mcdi - {$i/c)\\fj . 



9 

Proof . Consider the limit equation (|32p which can be written as 

(37) mc 2 F{x) - cG'{x) = XF(x) and cF'(x) - mc 2 G{x) = \G{x) . 

If i = j, then the result is obvious. So let i ^ j, we treat the case i < j, where the proof for 
i > j goes through mutatis mutandis by using forward difference approximations for derivatives. 
Assume i < j, also we prove the first assertion of the lemma, the proof of the second assertion 
is similar. Consider the second part of ([3?]) for Xj 

(38) cF'( Xj ) - mc 2 G( Xj ) = XG(xj) . 
Using backward difference approximations for derivatives we can write 

^ F( Xj ) - F(xj) _ fj-fj ^ 



(39) 

k=i+l 

Substituting ([39]) in ([38]) completes the proof. 



-J2 h * 

k=i+l 



We continue the proof of Theorem 1, consider the dominant parts with respect to c, so let 
c — > oo in ([35]) and ([36]) and simplify to get 

n 

( 40 ) (( ~ a k - Vk$k)X + {cQk - m 2 c 2 'r] k 'd k )T j + (mc 2 a k + mc 2 n k ^ k )^ fj + 

k=l 

n 

+ 1^2 \ { T 3Qk®k ~ TjT} k )X + {mc 2 Q k 'd k - mc 2 r) k )Tj + (m 2 c 3 a k i9 k + cn k )^j gj = 



k=l 



and 



n 

(41) [y~^ ( ( - TjTjk + TjQ k -d k )\ + (mc 2 r] k - mc 2 Q k tf k )Tj + ( - cr? fc - m 2 c 3 <r fe i? fc )) fj + 



fc=i 



+ [y^ f ( - %i?fe - crfe) A + (mVr/fctffc - c£»fc)Tj + ( - mc 2 n k -d k - mc 2 a k )^j gj = 0. 



fc=i 



To make the derivation simpler, the following notations are introduced 

a = Efe=i a fc = Efc=i( _cr fc ~ m^k), b = cbi- m 2 c 3 b 2 = YJk=ii c Sk ~ m 2 c 3 ?] k ^ k ), 
d = mc 2 di = YJk=i mc2 (°"fc + Vk'&k), e = Yl=i e k = J2k=i(&k^k ~ Vk), 

'Qi = Ylk=i mc 2 (g k -d k -r] k ), ui = ?n 2 c 3 ^i + cuj 2 = Efc = i(m- 2 c 3 o- fc t? fc + cn k ). 



mc 



By these notations, equations (140p and (]4ip can be written as 

aX + brj + d &TjX + qrj + uj \ f fj 
erj X — qTj — uj aX — brj — d J \ gj 



(42) 
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Since fj and gj are not identically zero for all j, then we expect 

(43) det ( aX + bT] +d eT J X + qT i +u ) =0i 

\ eTjX — qrj — gj aX — brj — d J 

where det(-) is the determinant of matrix. After simplifying, equation (|43p leads to 

(44) \±{ Tj ) = ± 1 



l {bTj+d) 2 - {qTj+ujf 



a 2 — e 2 Tj 

By |20j . the only accumulation point for the eigenvalue for the radial Coulomb-Dirac operator 
in the vicinity of x at infinity is mc 2 . So, we like to have 

|A+ — mc 2 \ = 

m 2 c 4 {a 2 -e 2 T 2 ) = (brj + d) 2 - {qrj + u) 2 

= {cb\Tj — •m l c'b<iTj + mc 2 di) 2 — (mc 2 qiTj + m 2 c 3 uji + c^) 2 - 

Letting m = 1, dividing both sides by c 6 , and taking the limit as c — > oo, we get 

(45) b 2 2 rf - u 2 = 0. 

Substituting back the values of 62 and oj\, the desired consequence is obtained for the fixed j as 

n n 

(46) Tj = j ^ Gk'&k j ^2 11k ' &k 

k=l k=l 

The above result can be generalized for arbitrary j as 

n n 

(47) Tj = \ Y, n :,■'>,■ / Yl ty&i 



i=l i=l 



and this ends the proof. 



The /ip-cloud functions depend strongly on the dilation parameter pj. As pj gets smaller, i.e., 
Pj max{hj, hj + \\ (= hj+x for exponentially distributed nodal points), as the shape functions 
of MMs in general become closer to the standard finite element functions, see Figure 3. In this 
case the FEM stability parameter might be applicable for MMs [17| 

FEM MMs . . h . 

Tj > Tj , as pj — > n J+ i. 

On the other hand, one should be careful about the invertibility of the matrix M, i.e., we can 
not approach pj = /ij+i which makes M singular. In Lemma 2, we derive the stability parameter 
t jem £ Qr ^Yiq computation of the eigenvalues of the radial Dirac operator, H K , using the FEM 
with quartic spline. The proof of the lemma is simple and uses the same technique as of the 
theorem above, thus we directly utilize the result of this theorem with minor modifications. In 
Table 7, the result of applying rJ EM for stabilizing the /ip-cloud method with pj = 1.1 • hj+\ is 
obtained, the approximation is good enough and the spuriosity seems to be eliminated. But a 
difficulty arises, that is, the end of the spectrum (the spectrum tail) behaves in a strange way, 
which may be regarded as spurious solutions. 
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hp-cloud shape functions for exponentially distributed nodal points hp-cloud shape functions for exponentially distributed nodal points 




domain O domain Q 



hp-cloud shape functions for exponentially distributed nodal points hp-cloud shape functions for exponentially distributed nodal points 




domain O domain O 



Figure 3. PU /ip-clouds with different dilation parameters: pj = 4- hj+i (up to 
the left), pj = 2 • hj + \ (up to the right), pj = 1.5 • hj + \ (below to the left), and 
Pj = 1.2 • hj + \ (below to the right). Quartic spline is used as a weight function. 



Lemma 2. The FEM stability parameter for the computation of the eigenvalues of the radial 
Dirac operator using quartic splines as a basis has the form 

(aq\ fem _ 3 h (hj+i - hj) 

J 17 [h j+ i + hj) 

Proof . Consider the general formula derived in Theorem 1 

n n 

(49) Tj = | n .ri } ji j Yl Htfji 

i=l i=\ 

where flji is defined by ([30]) . and o~ji and r)ji are respectively the entries of the matrices Mooo and 
Mioo- Note that in the FEM with quartic spline basis functions, Mooo an d Mioo are tridiagonal 
matrices with j th row elements as in Table 2. 



Table 2. The element integrals of the matrices Mqoo and M 



100- 



Index 

Matrix 


3-1 


3 


3 + 1 


j th row of Mooo 


3 h 


%( h j + hj+i) 


^hj + i 


j th row of Afioo 


17 
70 





17 
70 
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By Substituting the values of <jji and rjji from Table 2 in (|49p and using the definition of flji, 
we get the desired consequence. ■ 

5. Results and discussions 

Since the main goal of this work is applying the /ip-cloud method with the stability scheme, 
most of the discussion (all figures and tables except Table 7) provided here will be about the 
main stability parameter r,- in (I3ip given in Theorem 1. However, only Table 7 sheds some light 
on the FEM stability parameter given by Lemma 2. This discussion takes a form of comparison 
with the main stability parameter. 

For point nucleus, the relativistic formula is used to compare our results 

(50) A. 



1+ z2 « 2 



(n r -l+Vii 2 -Z 2 a' 2 ) 2 

where a is the fine structure constant which has, in atomic unit, the value 1/c, and the orbital 
level number n r takes the values 1, 2, To ease performing the comparison, the exact eigen- 
values \ nr ,K and the positive computed eigenvalues are shifted by —mc 2 . All computations are 
performed for the Hydrogen-like Ununoctium ion, where the atomic number and atomic weight 
for the Ununoctium element are 118 and 294 respectively. Consequently, and since the electron 
in the Hydrogen-like Ununoctium ion admits relatively large magnitude eigenvalues, for better 
measuring of the approximation accuracy, through out all computations we shall use the relative 
error. To treat the singularity of the pure Coulomb potential at x = 0, extended nucleus is as- 
sumed by modifying the potential to fit the finite nuclear size. The modified Coulomb potential 
considers another distribution of the charge along the nucleus (in the region [0 , R] where R is 
the nucleus radius) and pure Coulomb potential in the rest of domain. The continuity and the 
smoothness property (at least C 1 ) should be saved for the total modified potential. For the 
distribution of charge along nucleus, one can consider, e.g., uniform or Fermi distributions, in 
this work we consider uniformly distributed charge. 

As for the boundary conditions, the homogeneous Dirichlet condition is assumed. Note that 
for better approximation of the eigenstates lsi/2 and 2pi/ 2 , suitable Neumann boundary con- 
ditions should also be considered, see [2]. However, here, we do not treat these cases, instead, 
general computations are performed to account for the essence of discussion. The homogeneous 
Dirichlet boundary condition is then simply implemented, after coupling with the FEM, by 
omitting the two finite element functions at the lower and upper boundaries. 

As mentioned before, the computation of the radial Dirac operator eigenvalues requires ex- 
ponential distribution of the nodal points in order to capture desired behavior of the radial 
functions near the origin. For this purpose, we shall use the following formula 

(51) Xl = exp (ln(J a + e) + ( Hh + " ] ~ HIa + e) )*) - e , i = 0, 1, 2, . . . , n, 

where n is the total number of nodal points and e £ [0 , 1] is the nodes intensity parameter. The 
role of e is to control the intensity of the nodal points close to origin, as smaller e as more nodes 
are dragged toward the origin, see the discussion below. As for other approximation methods, 
increasing the number of nodal points provides better approximation, but this, of course, on the 
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account of the computational time. However, we can still obtain a good approximation with 
relatively less time, compared with increasing the nodal points, if the number of integration 
points is increased (the same size of the generalized matrices is obtained for a fixed number of 
nodal points, where increasing the number of integration points means more time is needed for 
functional evaluations but the same time is used for eigenvalues computation of the generalized 
system). This does not mean that we do not need to increase the number of nodal points to 
obtain more computed eigenvalues and to improve the approximation, but to get a better rate of 
convergence with less computational time, increasing both the numbers of integration and nodal 
points are necessary. In this computation, we fix the number of integration points at 10 ■ n. 

Table 3 shows the approximated eigenvalues of the electron in the Hydrogen-like Ununoctium 
ion obtained using the usual and the stabilized /ip-cloud methods, the computation is obtained 
at pj = 2.2hj+i, e = 10~ 5 , and n = 600. The clouds are enriched by P t {x) = [1, x(l — 
x/2) exp(— x/2)]. The eigenvalues, through out the computations in this work, are given in 
atomic unit. In Table 3, with the usual /ip-cloud method, the instilled spurious eigenvalues 
appear for both positive and negative k (the two shaded values in the fourteenth level), also 
the unphysical coincidence phenomenon occurs for the positive k (the shaded value in the first 
level). Note that these spuriosity of both categories are removed by the hp-CP G method. 

Table 3. The first computed eigenvalues of the electron in the Hydrogen-like 
Ununoctium ion using the usual and the stabilized /ip-cloud methods for point 
nucleus. 



Level 


/ip-cloud 


ftp-cloud 


Exact solution 


ftp-CPG 


ftp-CPG 




K = 2 


n = -2 


K = -2 


K = -2 


K = 2 


1 


-1829.630750899 


-1829.630750902 


-1829.630750908 


-1829.628309112 




2 


-820.7698136330 


-826.7698136329 


-826.7683539069 


-826.7714785272 


-826.7738882959 


3 


-463.1214970564 


-463.1214970566 


-463.1183252634 


-463.1247150569 


-463.1261170024 


4 


-294.4552367950 


-294.4552367952 


-294.4509801141 


-294.4591541031 


-294.4600671778 


5 


-203.2468937049 


-203.2468937047 


-203.2419549027 


-203.2511517040 


-203.2517946674 


6 


-148.5588260984 


-148.5588260983 


-148.5534402360 


-148.5632453116 


-148.5637243357 


7 


-113.2536099083 


-113.2536099084 


-113.2479180697 


-113.2580871797 


-113.2584595495 


8 


-89.16385480233 


-89.16385480237 


-89.15794547564 


-89.16832365853 


-89.16862284813 


9 


-72.00453396071 


-72.00453396065 


-71.99846504808 


-72.00894720487 


-72.00919403005 


10 


-59.35481340095 


-59.35481340100 


-59.34862423729 


-59.35913470352 


-59.35934276227 


11 


-49.76429096817 


-49.76429096819 


-49.75800915710 


-49.76849047005 


-49.76866900765 


12 


-42.32147184311 


-42.32147184312 


-42.31511730902 


-42.32552373918 


-42.32567925216 


13 


-36.43039621976 


-36.43039621984 


-36.42398370073 


-36.43427738957 


-36.43441456989 


14 


-33.96502895994 


-33.96502895893 


-31.68173025393 


-31.69187884728 


-31.69200116063 


15 


-31.68818961940 


-31.68818961935 


-27.80813459180 


-27.81810976712 


-27.81821982418 



5.1. Integration of /ip-cloud functions 

To approximate the integrals in the weak form in the Galerkin hp-cloud approximation, we 
use two-point Gaussian quadrature rule. Gaussian quadrature rules are the most used numeri- 
cal techniques to evaluate the integrals in MMs due to their exactness property in integrating 
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of polynomials of degree 2m q — 1, where m q is the number of quadrature points [18]. However, 
using Gaussian quadrature rules yields integration error when the grids are not coincident with 
the clouds covers, and thus instabilities and spurious modes start to appear. Also for non- 
uniformly distributed points (the case we assume in this work), Gaussian quadrature rules do 
not pass the patch test (fail in consistency). Therefore, stabilizing conforming nodal integration 
(SCNI), see [9], is introduced to overcome these difficulties. The main feature of SCNI is using 
the divergence theorem to substitute the derivative, i.e., the derivative ^ i ^ h in the sub-domain 
Vtj = [xj , xj+i] is replaced by a smooth derivative (averaging derivative) 4z:^ h at x € £lj as 



-*»(z) » ~^{x) = -rmx)dx 



dx dx Xj-^-~^ ^ j J x ■ dx ^ j 

This definition helps stabilizing the integration, further, it saves time in the computation by not 
calculating the derivatives of the cloud functions. Thus, there is no need to evaluate (M -1 )' = 
—M M'M -1 , which is expensive to calculate. For integrating and programming the weak form 
in MMs, the results from [131 Q3] are useful. 

The cloud shape functions are evaluated at the integration points (digital evaluation), since 
, practically, it is somehow impossible to write the cloud functions explicitly without matrix 
inversion. Also, it is not recommended to obtain the inverse of M directly, instead, LU factor- 
ization is better to be used from cost (less time consumption) and numerical stability point of 
views. Moreover, in MMs generally, to enhance the stability of the computation and to main- 
tain the accuracy (that may be affected or lost due to the round-off error), and to get better 
conditioning of the matrix M (lower condition number), the origin should be shifted to the 
evaluation point |18[ 1241 127] . i.e., x is replaced by the transformation x — x Xqy%q} conse 

n 

quently ipi(x) = P t (0)M~ 1 (x)B i (x) where M(x) = y~]ipj( -)P(xj- x or i g )P t (x i - x orig ) and 

i=i Pi 

Bi(x) = i Pl (^)P(x i - X ori g). 

5.2. Enrichment basis functions P(x) 

For the reason discussed before, only intrinsic enrichment, P(x), is considered in the definition 
of the /ip-cloud functions for the computation of the eigenvalues of the radial Dirac operator. 
The number and the type of enrichment functions in the basis set P(x) can be chosen arbitrary 
for each cloud \19\ [32] , but for practical reasons (lowering both the condition number of M and 
the computational time) we assume P(x) = [l,pi(x)\. For the approximation of the radial Dirac 
operator eigenvalues, to enrich the cloud with a suitable basis P(x), two main properties should 
be considered; firstly, and sufficient one, the elements of P(x) ought to have the continuity prop- 
erties (continuous with continuous first derivatives) of the space "K so that for all j, the cloud ipj 
is a C 1 -function, provided that (fj is a C 1 -function. Secondly, global behavior and fundamental 
characters about the electron motion of the Hydrogen-like ion systems should be embedded in 
P{x). Slater type orbital functions (STOs) and Gaussian type orbital functions (GTOs) pro- 
vide good description of the electron motion |10[ [21] . The quadratic term in the exponent of 
the GTOs causes numerical difficulty, that is, with the GTOs, the matrix M rapidly becomes 
poorly conditioned, this is also what is observed when applying quadratic basis enrichments, see 
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[6]. Consequently, the STOs are considered as the enrichment of the /ip-cloud functions, thus 
Pi(x) can have, e.g., the following forms 

exp(— x), xexp(—x/2), x{l — x/2) exp(— x/2), . . . etc. 

Note that, these functions should be multiplied by normalization parameters, but, computation- 
ally, there is no effect of multiplication by these parameters. 

Since the global behavior of the eigenstates of the Hydrogen-like ions in the relativistic case 
(Dirac operator) does not differ much from that of the non-relativistic case (Schrodinger oper- 
ator), one can also assume the solutions of the radial Coulomb-Schrodinger eigenvalue problem 
as intrinsic enrichments (see e.g. [22]) 

%n r i( x ) = Nn r £ (2Zx/n r a o y L^^ e (2Zx/n r a ) exp(-Zx/n r a ), 

where h 2e+ .\(x) = - — — f Ur ^ ^ ] x k [ s the Laguerre polynomial, an is the Bohr 
- k - V n r +£- k J & ^ J 

radius, n r = 1,2,... is the orbital level number, and I is, as mentioned before, the orbital 

angular momentum number given to be zero for s-states, one for p-states, two for d-states, etc. 

For general intrinsic enrichment, it is, somehow, tedious to apply the above formula for each 

level n r , instead, good results are still achievable even with, e.g., n r equals the first possible level 

of the given state (i.e., n r = 1 for all s-states, n r = 2 for all p-states, n r = 3 for all o!-states, 

etc.). Moreover, it is also possible to consider enrichment based on the solution of the radial 

Dirac eigenvalue problem, see e.g. but the above enrichments are simpler from practical 

point of view. In the coming discussion, the enrichment basis P t {x) = [1 , x(l — x/2) exp(— x/2)] 

is assumed in all computations. 

5.3. Dilation Parameter p 



The dilation parameter, p, plays a crucial role in the approximation accuracy and stability, 
it serves as the element size in the FEM. The parameter p can be chosen fixed or arbitrary, but 
it is often assumed to be constant for all /ip-clouds when uniformly distributed nodal points are 
used. In this work, exponentially distributed nodal points are assumed to get enough informa- 
tion about the radial functions near the origin where they oscillate heavily relative to a region 
away from it. Thus we consider 

Pj = v ■ max{hj , /ij+i} = uhj + i, 

where the maximum is considered to engage sufficiently large region where the cloud function 
is defined so that the risk for singularity of the matrix M is substantially decreased, further, v 
is the dimensionless size of influence domain [27] . Moreover, for a non-uniform mesh, v can be 
chosen locally, i.e., v = vj, where, in this work, we assume fixed v. Now it remains to determine 
the value/values of v taking into account that pj should be large enough (y > 1) in order to 
ensure the invertibility of M (to ensure that any region is covered by at least two clouds). On 
the other hand, pj should not be very large to maintain local character of the approximation. As 
discussed before (see also Figure 3, the case v = 1.2), if v — > 1, then ijjj will act as finite element 
shape function, and thus the features of the /ip-clouds are gradually lost, also large values of v 
make tjjj to behave like polynomial of higher degree (see Figure 3, the case v = 4). To conclude, 
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v should be chosen moderately and such that it guarantees that no integration point is covered 
by only one cloud [271 E2] • 

The optimal choices of v can be determined individually for each problem by carrying out 
numerical experiments. In [30^ 144] . it is shown that v € [2,3] provides nice results for the elas- 
ticity problem. For the computation of the radial Dirac operator eigenvalues with the stability 
scheme, when v € [2.2, 2.7] good results are obtained and the spurious eigenvalues are completely 
eliminated. Also as v gets smaller as better approximation is obtained, see Table 4. In Figure 

Table 4. The first computed eigenvalues of the electron in the Hydrogen-like 
Ununoctium ion for k = —2 for point nucleus with different values of u, where 
n = 600 and e = 10 -5 are used. 



Level 


v = 2.0 


v = 2.2 


v = 2.5 


v = 2.7 


Exact solution 


1 


-1829.6287 


-1829.6283 


-1829.6276 


-1829.6270 


-1829.6307 


2 


-826.77119 


-826.77147 


-826.77197 


-826.77233 


-826.76835 


3 


-463.12417 


-463.12471 


-463.12567 


-463.12638 


-463.11832 


4 


-294.45850 


-294.45915 


-294.46033 


-294.46120 


-294.45098 


5 


-203.25046 


-203.25115 


-203.25244 


-203.25340 


-203.24195 


6 


-148.56255 


-148.56324 


-148.56460 


-148.56562 


-148.55344 


7 


-113.25741 


-113.25808 


-113.25949 


-113.26054 


-113.24791 


8 


-89.167688 


-89.168323 


-89.169756 


-89.170831 


-89.157945 


9 


-72.008358 


-72.008947 


-72.010396 


-72.011489 


-71.998465 


10 


-59.358602 


-59.359134 


-59.360592 


-59.361700 


-59.348624 


11 


-49.768025 


-49.768490 


-49.769950 


-49.771070 


-49.758009 


12 


-42.325133 


-42.325523 


-42.326981 


-42.328113 


-42.315117 


13 


-36.433970 


-36.434277 


-36.435728 


-36.436870 


-36.423983 


14 


-31.691663 


-31.691878 


-31.693318 


-31.694472 


-31.681730 


15 


-27.817992 


-27.818109 


-27.819533 


-27.820699 


-27.808134 



4, we study the convergence rate of the first five eigenvalues in Table 4. 

The effect of the influence domain factor, v, on the convergence 

-4.2 



-4.6' 



o -5.4 
-5.6 
-5.8 
-6' 

2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 

influence domain factor v 

Figure 4. Studying the convergence rate with respect to the influence domain 
factor v. The comparison is carried out for the first five eigenvalues in Table 4. 

It is clear how the smaller v gives the better approximation. One argues, as it is clear from 
the figure, that v can be, e.g., of some value less than 2 in order to achieve a better rate of 
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convergence. However, this will be on the account of spuriosity elimination (the cloud is not 
stretched enough to capture the desired behavior of the approximated solution) and on the ac- 
count of the invertibility of the matrix M (for small v some regions are covered with one cloud) . 
However, as in the FEM, one can apply /i-refinement in the /ip-cloud method (see e.g. [161 144| ). 
this can be done by assuming smaller values of the dilation parameter pj (keeping v fixed and 
making hj + \ smaller by increasing the number of nodal points). Thus as pj getting smaller, 
more clouds of smaller domain sizes are added. 

The intensity of the exponentially distributed nodal points near the origin has an influence 
on the convergence rate of the approximation. The intensity of the nodes, near the origin or 
away from it, is controlled by the nodes intensity parameter, e, via formula (|5ip . As smaller 
value of e is considered as more concentration of nodes near the origin is obtained, see Figure 5 
(the graph to the left). 

Table 5 shows the computation of the eigenvalues with different values of e with 600 nodal 
points. The computation with e smaller than 10~ 7 is almost the same as of e = 10~ 7 , thus it is 
not required to study smaller values of e than 10 -7 . 

Table 5. The first computed eigenvalues of the electron in the Hydrogen-like 
Ununoctium ion for k = — 2 for point nucleus with different values of e, where 
n = 600 and v = 2.2 are used. 



Level 


e = l(r 4 


e = 10~ 5 


e = 10~ 6 


e = 10- 7 


Exact solution 


1 


-1829.6289 


-1829.6283 


-1829.6280 


-1829.6280 


-1829.6307 


2 


-826.77073 


-826.77147 


-826.77170 


-826.77173 


-826.76835 


3 


-463.12322 


-463.12471 


-463.12517 


-463.12523 


-463.11832 


4 


-294.45726 


-294.45915 


-294.45973 


-294.45981 


-294.45098 


5 


-203.24904 


-203.25115 


-203.25180 


-203.25188 


-203.24195 


6 


-148.56101 


-148.56324 


-148.56393 


-148.56402 


-148.55344 


7 


-113.25578 


-113.25808 


-113.25879 


-113.25888 


-113.24791 


8 


-89.165992 


-89.168323 


-89.169039 


-89.169131 


-89.157945 


9 


-72.006610 


-72.008947 


-72.009662 


-72.009755 


-71.998465 


10 


-59.356811 


-59.359134 


-59.359844 


-59.359936 


-59.348624 


11 


-49.766195 


-49.768490 


-49.769189 


-49.769279 


-49.758009 


12 


-42.323268 


-42.325523 


-42.326208 


-42.326296 


-42.315117 


13 


-36.432073 


-36.434277 


-36.434943 


-36.435030 


-36.423983 


14 


-31.689734 


-31.691878 


-31.692524 


-31.692607 


-31.681730 


15 


-27.816033 


-27.818109 


-27.818732 


-27.818812 


-27.808134 



In Figure 5 (the graph to the right), the first computed eigenvalues of Table 5 are studied. 
It is clear that as e gets larger (up to some limit), better approximation is obtained. However, 
as mentioned before, the rate of convergence is almost the same when e € (0 , 10~ 7 ) (e = is 
excluded to avoid ln(0) when extended nucleus is assumed). Also e does not admit relatively 
large values in order to get enough nodes close to the origin, where the radial functions oscillate 
relatively more, without increasing the number of nodal points. Therefore, the most appropriate 
values for e which provide good results, are somewhere in the interval [10 -6 , 10 -4 ]. 
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Different exponential distributions ol fifty nodes The effect of file exponential distributions on the convergence 




Figure 5. To the left, different exponentially distributed nodal points are plot- 
ted using the formula (|5ip . To the right, the effect of nodes intensity near the 
origin on the convergence rate, the comparison is carried out for the first five 
eigenvalues in Table 5. 

The approximation of the stabilized /ip-cloud scheme with different numbers of nodal points 
is given in Table 6. The rate of convergence of the corresponding first five eigenvalues is studied 
in Figure 6, where h is the maximum of all distances between the adjacent nodes which equals 
to h n = x n — x n -i, the distance between the last two nodes for exponentially distributed nodes. 

Table 6. The first computed eigenvalues of the electron in the Hydrogen-like 
Ummoctium ion for k = —2 for point nucleus with different number of nodes, 
where v = 2.2 and e = 10~ 5 are used. 



Level 


n = 200 


n = 400 


n = 600 


n = 800 


n = 1000 


Exact solution 


1 


-1829.5628 


-1829.6224 


-1829.6283 


-1829.6297 


-1829.6302 


-1829.6307 


2 


-826.82670 


-826.77726 


-826.77147 


-826.76987 


-826.76923 


-826.76835 


3 


-463.23292 


-463.13630 


-463.12471 


-463.12146 


-463.12016 


-463.11832 


4 


-294.59147 


-294.47367 


-294.45915 


-294.45503 


-294.45336 


-294.45098 


5 


-203.39386 


-203.26721 


-203.25115 


-203.24654 


-203.24466 


-203.24195 


6 


-148.70878 


-148.58009 


-148.56324 


-148.55835 


-148.55635 


-148.55344 


7 


-113.40170 


-113.27527 


-113.25808 


-113.25304 


-113.25096 


-113.24791 


8 


-89.306709 


-89.185557 


-89.168323 


-89.163201 


-89.161076 


-89.157945 


9 


-72.139617 


-72.026008 


-72.008947 


-72.003802 


-72.001653 


-71.998465 


10 


-59.480154 


-59.375861 


-59.359134 


-59.354006 


-59.351849 


-59.348624 


11 


-49.878353 


-49.784751 


-49.768490 


-49.763410 


-49.761256 


-49.758009 


12 


-42.423104 


-42.341207 


-42.325523 


-42.320517 


-42.318374 


-42.315117 


13 


-36.518814 


-36.449288 


-36.434277 


-36.429365 


-36.427242 


-36.423983 


14 


-31.762955 


-31.706134 


-31.691878 


-31.687081 


-31.684984 


-31.681730 


15 


-27.875610 


-27.831538 


-27.818109 


-27.813442 


-27.811376 


-27.808134 



The lack of error estimates for the approximation of the Dirac eigenvalue problem due to the 
boundedness problem results in an incomplete picture about the convergence analysis. Never- 
theless, from Figure 6, the convergence rate of the approximation of the first five eigenvalues, Ai, 
. . ., A5, is nearly 3.09, 2.66, 2.62, 2.59, and 2.56 respectively, which it takes a slight decreasing 
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pattern as we go higher in the spectrum levels, see the corresponding table. 




Figure 6. Study the convergence rate of the first computed five eigenvalues in 
Table 6. 

With the stability parameter r FEM , the computation is presented in Table 7. The computation 
is obtained with 600 nodal points at v = 1.1 and e = 10~ 5 . The result is compared with the 
same stability scheme but with the stability parameter r at the same parameters but v = 2.2, 
the comparison is also obtained in the non-relativistic limit (very large c). 

Table 7. The first computed eigenvalues of the electron in the Hydrogen-like 
Ununoctium ion for k = — 2 for point nucleus using the stability scheme with the 
stability parameters r and t fem . 







The speed of light 






100 x The speed of light 




Level 


T 


T FEM 


Exact solution 


r 


T FEM 


Exact solution 


1 


-1829.6283 


-1829.6304 


-1829.6307 


-1740.2372 


-1740.4777 


-1740.5080 


2 


-826.77147 


-826.76993 


-826.76835 


-773.73860 


-773.57259 


-773.56033 


3 


-463.12471 


-463.12174 


-463.11832 


-435.46054 


-435.14787 


-435.12752 


4 


-294.45915 


-294.45551 


-294.45098 


-278.88245 


-278.49775 


-278.48144 


5 


-203.25115 


-203.24715 


-203.24195 


-193.82362 


-193.39522 


-193.38978 


6 


-148.56324 


-148.55905 


-148.55344 


-142.53145 


-142.07261 


-142.08222 


7 


-113.25808 


-113.25377 


-113.24791 


-109.23625 


-108.75140 


-108.78165 


8 


-89.168323 


-89.163916 


-89.157945 


-86.404375 


-85.894382 


-85.950912 


9 


-72.008947 


-72.004478 


-71.998465 


-70.067886 


-69.534118 


-69.620219 


10 


-59.359134 


-59.354644 


-59.348624 


-57.975599 


-57.420335 


-57.537357 


11 


-49.768490 


-49.764010 


-49.758009 


-48.773149 


-48.197640 


-48.347352 


12 


-42.325523 


-42.321064 


-42.315117 


-41.606088 


-41.009232 


-41.195370 


13 


-36.434277 


-36.429826 


-36.423983 


-35.913753 


-35.292096 


-35.520492 


14 


-31.691878 


-31.687405 


-31.681730 


-31.315908 


-30.664671 


-30.942291 


15 


-27.818109 


-27.813579 


-27.808134 


-27.547311 


-26.861865 


-27.195369 



As it can be seen from Table 7, the convergence property with 

T FEM ig 

slightly better. Unfor- 
tunately, the approximation with t fem seems to behave strangely at the end of the spectrum, 
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that is, only the spectrum tail has the following behavior (the last eigenvalues of the computa- 
tion in Table 7 with r M for the relativistic case) 



A- 



mc 



A_ + mc 2 



207072481.0215 
211429663.4158* 



-215565247.3448 
-220006205. 1800* 



226003907.3130 
231896256.0483* 



-235294474.7992 
-241138935.9851* 



246890583.9362 

257292411.7094* 

267659710.2673* 



-257366374.4374 

-267386241.2969* 

-279193268.7275* 



291928112.6166 
296228215.8873* 



-303237209.5231 
-308029351.9019* 



This behavior occurs only for few values at the end of the spectrum, and no such effect is revealed 
in the rest of the spectrum. To our knowledge, the values marked with * might be spurious 
eigenvalues for some unknown origins in higher levels, which, in calculating the correlation 
energy, seem to have no significant effect. 

Table 8 shows the computation of the eigenvalues of the electron in the Hydrogen-like Un- 
unoctium ion with k = —2. The computation is for extended nucleus obtained using the stability 
scheme, where the first and the last computed eigenvalues are presented. The number of nodes 
used is 1000, also the used values of u and e are respectively 2.2 and 10 -5 . 

Conclusion. 

The scheme developed in this work, the hp-CP G method, for stabilizing the /ip-cloud ap- 
proximation for solving the single-electron Coulomb-Dirac eigenvalue problem ensures complete 
treatment of the spurious eigenvalues. The scheme strongly depends on the derived stability 
parameter r, which is simple to implement and applicable for general finite basis functions. 
The elimination of the spurious eigenvalues is also affected by the influence domain factor v, 
for v less than 2, spurious eigenvalues start to appear. The convergence rate is high for the 
first eigenvalues, while it slowly decreases as the level gets higher. Comparing with the finite 
element stability approach [2], the scheme convergence rate is lower. We may state that, as the 
main disadvantage of MMs in general, the /ip-cloud method is more expensive due to the time 
consumption in evaluating the shape function which demands more integration point as v gets 
larger to obtain the desired accuracy. The number of integration points used here is ten times 
the number of nodal points (this large number of points is assumed in order to study the effects 
of the other parameters from a comparative point of view), which can be made smaller, i.e., 
v > 2 is enough to get sufficient accuracy. 
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Table 8. The first and the last computed eigenvalues of the electron in the 
Hydrogen-like Ununoctium ion for k = —2 for extended nucleus using the stability 
scheme. 
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